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Abstract 

The recent interest in the imaging possibihties of photonic crystals (superlensing, superprism, 
optical mirages etc..) call for a detailed analysis of beam propagation inside a finite periodic 
structure. In this paper, we give such a theoretical and numerical analysis of beam propagation 
in ID and 2D photonic crystals. We show that, contrarily to common knowledge, it is not always 
true that the direction of propagation of a beam is given by the normal to the dispersion curve. We 
explain this phenomenon in terms of evanescent waves and we construct a renormalized dispersion 
curve that gives the correct direction. 

PACS numbers: 42.70 Qs, 42.25.Fx 



I. INTRODUCTION AND SETTING OF THE PROBLEM 



Some beautiful experiments and numerical works have shown that it was possible to 



obtain quite unusual behaviors of 
crystals (PhCs). 



ight propagation inside meta-materials and photonic 



]. In particular, recent ideas by Pendry confirmed by 
experiments show that photonic crystals, maybe under the guise of meta-materials, could 
prove to be of huge interest in that they could allow to beat the diffraction limit and to make 
superlenses. The point of this work is to give a theoretical insight into beam propagation 
inside PhCs and in particular on the computation of the shift of the transmitted beam (cf. 
Figure 1). We consider both ID and 2D PhCs and we show that, contrarily to what is 
generally believed, the direction of propagation is not always directly given by the normal to 
the dispersion curves for 2D PhCs. Rather, we define a renormalized, or effective, dispersion 
diagram, whose normal gives the correct direction of propagation. Numerical examples are 
given illustrating the various regimes. 

We will begin by studying, systematically, ID structures before extending the results to 
2D structures (seen as stacks of gratings) through the concept of the so-called two waves 
approximation that will be introduced later. Throughout this work, we use time-harmonic 
fields, with a time-dependence of exp(— i0t), that are ^-independent. The vectorial diffrac- 
tion problem is reduced to the study of the two usual cases of polarization: s-polarization 
(electric field parallel to the grooves of the gratings) or p-polarization (magnetic field parallel 
to the grooves). The wavenumber is denoted by /cq = X' where A is the incident wavelength 
in vacuum. 

The incident field is a Gaussian beam whose z component can be expressed by 

u' (x, y) = ^ A (a) e<""+^/^^)rfa (1) 

where 

A(«) = ^e-^(— 0^ 

we denote the electric (magnetic) field in the case of s (p) -polarization and ao = ko sin 6*0, 
where 6o is the mean angle of incidence of the beam and w its waist. 
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II. ANALYSIS OF THE BEAM PROPAGATION 



A. The case of a stratified medium : ID PhC 

In this section, we derive the value of the shift of the transmitted beam for the particular 
case of a stratified medium, i.e. when the relative permittivity is constant in the horizontal 
direction: it is described by a real periodic function (period h): e{y). We denote: r = 
{x,y). We consider N periods of the stratified medium, which is embedded in vacuum. 
For an incident plane wave of wavevector k = A;(sin^, — cos^^,0), we denote /?o = ^ocos^ 
and {rfyf {k,6) ,t]sf {k,6)) the reflection and transmission coefficients of the structure, the 
electromagnetic field in the outer regions reads as: 

u{r) = e*''-" + r^e''''''^''''''^+y^°'^\ y>0 (2) 
u{r) = t^e^'^o(^^''^^-(^+^'^)^°'^), y < -Nh (3) 

We denote by T the transfer matrix of one period, then it is known 91] that the refiection 
and transmission coefficients are related through the relation: 

Let us denote by 7 and 7"^ the eigenvalues of T and by v = (^n, 02i) , w = (0i2, 022) the 
associated eigenvectors (Tv = 7V , Tw = 7^^w). It is known (see for instance 91]) that the 
band gaps and the conduction bands are given respectively by: G = {{k,6) ,tr (T) > 2}, 
and : B = {{k,9) ,tr (T) < 2}. The reflection and transmission coefficients are then given 
by: 

^) = !yL _ gif ' ^) = (5) 

where, denoting q{x,y) = ^^j^, the functions / and g are deflned by g{k,6) = 
q{\'),f{k,6) = g (w) and v is chosen such that \g\ < 1 in the conduction bands. Re- 
mark that in these bands, the inverse of / is equal to the conjugate of g (see 10|, lUl for 
details). 
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A straightforward calculation shows that, for {k, 9) e B: 



+ (X} 



T^iKO) = ^ + (6) 

+00 

p=0 

In the conduction bands, the eigenvalue 7 can be written under the form: 7 = e^^^, where 
(3 is the so-called Bloch phase. When the incident field is a beam, we get an infinite sum 
of transmitted beams, corresponding to multiple scattering. Let us concentrate on the first 
transmitted beam, i.e. the beam that reads: 



,Nh)^ J A (a) (1 - 1^1^) e^^'^^e^^^da (8) 



whose Fourier transform is: 

Ut{a)^V2^A{a)e'>^'"' (9) 

where 

A{a) = A{a){l-\gf). (10) 

We denote Gi,Gt, Gd the points where, respectively, the incident, transmitted and re- 
flected beams enter or emerge from the PhC. Given the incident field, the axis are chosen 
to have Gi = 0. These points are defined as the barycenters, or first moments, of the 
corresponding fields, that is: 

„ / x M*(a::,0) dx 

" ^\ud{x,Nh)\^dx ^ ' 



G, 



j x\u (.r./i/i) ax 
1^. 



f\u*{x,Nh)\^dx 

Using Parseval-Plancherel identity, we get the angular shift due to the beam propagation 
(cf. fig. 1): 

Gt [ [a] dad (a) da , , 

J A' (a) da ^ ^ 

A series expansion of tan ip can be obtained provided the phase function is analytic with 
respect to a in a neighborhood of ckq. Indeed, we can then write: 

daP (a) = V ^^^liM. _ (13) 
^-^ ml 

m=0 
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We obtain after some manipulations: 

tanV^ = -^ '"'^^^;^+^^ gr-^VM (14) 

where F is the Euler Gamma function [12]. When kow is large, then A {a) is concentrated 
around Oq, and if daP (a) does not vary too quickly in the vicinity of ao, we obtain the 
well-know crude approximation: 

tanip ~ —daP («o) (15) 

Of course, the formula (fT5|) can no longer hold if da(3 (a) is not analytic near ao, i-e. when 
ao is a branch point. We shall encounter this case in the following section. 

In order to give a geometric interpretation of this result, let us remark that {daP (ao) , —1) 
is a vector that is normal to the dispersion curve at wavelength A. So that we retrieve the 
well-known fact that for a spatially large beam, the direction of propagation is given by the 
normal to the isofrequency Bloch diagram. 

We shall see in the following that this result is in general no longer true in finite 2D 
structures. 



B. Beam propagation in a 2D photonic crystal 

The crystal is described as a stack of gratings and we assume that in the spectral domain 
defined by the above beam, the ratio between the wavelength and the period d of the gratings 
is such that there is only one reflected and one transmitted order for the grating structure. 
Then the propagating reflected and transmitted fields can be expressed as: 

(x, y) = [ A (a) (a) e'^'^^-'^yUa (16) 



ut {x, y)= A (a) t^ (a) e'^^'+'^^^rfa (17) 



Once the reflection and transmission coefficients (r7v,t7v) are known (by using a rigorous 



numerical method, for instance the Fourier Modal Method (FMM) [13|) there exists a unique 
unimodular real 2x2 matrix Tjy [9] with real coefficients satisfying relation (jl]): it is the 
dressed transfer matrix of the total structure We have: 



■TV 




-1 




e-'^^f^^ I \ 021 <p22 



(18) 



where the phase /?7v is the renormahzed Bloch phase for the global structure pi|. From this 
matrix, the reflection and transmission coeflicients can be written in the following form: 

rN = ^ tN= ~ ^ ^ ^' (19) 

Q2i^NNh _ j Q2i^NNh _ g-1 J ^ ' 

For a stratifled medium with homogeneous layers, the reduced transfert matrix satisfles 
rigorously the relation: = Tjv for all A^. However, for a two dimensional photonic 
crystal, this relation tends to become false as the number of periods is increased: this is due 
to the fact that matrix Ti does not take the evanescent waves into account. Consequently, 
as the thickness of the device increases the discrepancy between and Tjv increases as 
well. This remark has a crucial importance for our study as, in general, the derivative of 
the phase da^N is not equal to daP. 

By deflnition tr (Tat) = 2 cos (^Nhp^^, so that: 



djr (Tjv) = -2Nhdo,l3N sin [nHPn) = TNhd^(3N\/^ ~ tr (T 
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N) 



N 



the sign is unambiguously 



this provides us with a numerical method for computing 
flxed using the fact that < 1 and is associated with the eigenvalue e^^^P^ _ 

We have reduced the problem of computing the transmitted fleld to the one-dimensional 
case, and thus we can write: 



u'^{x,Nh)= A{a)e'^'"^^^°'^e"'''da (20) 



We can now give the main result of this paper, whose proof is given in Appendix 3. 
We assume that the beam is spatially large (i.e. kow ^ 1). Then two cases may be 
encountered with respect to the dispersion curve. For the mean angle of incidence of the 
beam (corresponding to ao), the curve is either regular (i.e the slope is not infinite ), or it 
is singular, i.e. the tangent to the curve is vertical. The shift of the beam is then described 
accordingly: 



1. If the tangent is not vertical, i.e. 



< +00 then the angle of refraction ijj of 



the beam inside the structure is given by: 

tamjj dapNiao), (21) 
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2. if the slope is infinite, i.e. 



+00 then the angle of refraction ip of the 



beam inside the structure is given by: 



23/4 /g\ 

-C^-^r (^-j Vsin^i + sin^oV^, (22) 

where ol\ = fcosin6'i is the maximum of A (a)^ {a — ao)^''^ and Cat is a constant such 
that /^AT ~ Cjq\/ — cto near ao- 

In the second result, the geometry of the structure and its electromagnetic parameters 
enter in the constant Cat. For a sufficiently large beam, Q\ ^ Oq. 

Two important properties should be noted in that case. The obvious one is that the shift 
does not tend to infinity when the normal to the dispersion curve tends to the horizontal 
axis, a fact that was of course expected, but which shows that the normal to the renormalized 
dispersion curve gives the direction of the beam, only if d^P does not vary too quickly in the 
vicinity of the mean angle aQ. Two parameters are in fact needed for a complete description 
of the situation: the normalized waist kow and the derivative of the phase daP (ao). The 
above result only gives the asymptotic behavior for the separate parameters. 

The second important point is the dependence of the shift with respect to the normalized 
waist, a situation which was not encountered in the first case. In order to understand this 
point, one should note that there is here a guided mode in the structure, i.e. a pseudoperiodic 
mode whose wavevector has a null vertical component. Of course, for a finite size structure, 
the uniqueness of the scattering problem implies that guided modes are associated with 
complex values of a which are poles of the transmission coefficient. The finiteness of the 
structure provokes a splitting of the eigenvalue into two complex values jsl corresponding 
to a zero and a pole of the refiexion coefficient. When such a structure is illuminated by 
a plane wave under the incidence ao the transmission shows a Fano profile indicating the 
excitation of the lossy mode. When the incident light is a beam, the behavior of the field 
resembles that of a plane wave in the limit k^w ^ 1, therefore the displacement of the 
barycenter towards infinity is associated with a spreading of the transmitted beam and 
thus, precisely because of the spreading, the very notion of barycenter of the transmitted 
beam loses its physical meaning. 
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III. NUMERICAL EXAMPLES 



In the following, we present some numerical computations illustrating the various situa- 
tions described by the above results. We will denote by: 

• A the shift computed by direct numerical computations of the fields. 

• the shift computed through the isofrequency dispersion diagram. 

• A^ the shift computed by use of the effective theory developed in the previous section. 

A. Case of a stratified one dimensional medium 

In this subsection, we check the numerical method that allows to compute the derivative 
of the Bloch vector of the equivalent T matrix and also the formula that gives the shift 
of the beam. The structure that we use is just a Bragg Mirror with two alternating slabs 
(thicknesses hi and /i2) in each period. The s polarized incident monochromatic beam is 
characterized by its waist w — 15A and its mean angle of incidence = 50°. The wavelength 
is such that \/hi = 2.27 with hi/h2 = 2 and the following permittivities for the slabs: 
El = 2.1, 82 = 6.25. In fig. 2, we give the amplitude of the incident, transmitted and reflected 
fields on the upper and lower interfaces of the device for = 15 periods. The shift of the first 
transmitted beam obtained by direct numerical computation is A. /hi — 12.51 whereas the 
shift obtained by computing the Bloch coefficient is Ap/hi — 12.52 and finally, the numerical 
computation described in the above section gives A^//ii — 12.52 hence, as expected, a 
perfect agreement with the Bloch approach. In order to complete this verification, we now 
use a p-polarized incident beam with w = 15A, X/hi = 4.75, = 47.5° and the parameters 
£i = 10.89, £2 = 1,^1/^2 = 1- Fig-3, shows the amplitude of the transmitted and reflected 
fields. This time, we obtain A/ hi = 62.7, A^//ii = 62.7209, A^//ii = 62.7209, which 
confirms the validity of our approach for that straightforward situation. 

B. Stack of gratings 

We are now in a position to apply our theoretical approach to the more complex situation 
of a stack of gratings. We recall that we assume that the wavelength is such that there is 
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only one Bloch mode inside the structure and one transmitted order and one reflected order. 
For all the following numerical experiments the field is s-polarized. 

The photonic crystal is a stack of 7 lamellar gratings with inverted contrast: Si = Sext = 
11.56 and 62 = 1( di/d = hi/d = 1/2 and h = d, see fig. 2 for notations). We compute 
the field for an incident beam {w = 12.5X,X/d = 2.2,6 = 50°) using the Fourier Modal 



Method 



13| . The amplitudes of the field on the upper and lower interfaces are given in 



fig.4. The Bloch diagram is given in fig.5. The shift of the transmitted beam obtained 
directly from this computation (through the envelope) is A/c? ~ 11.25 the shift computed 
from the isofrequency dispersion diagram is Ap/d = 2.45 and the shift obtained from the 
effective theory is A^/rf = 11.4. Therefore, we see that we have an error factor of 4.5 by 
neglecting the evanescent waves. 

The effective theory also applies when contra-propagative Bloch modes exist in the struc- 
ture, these modes authorizing super-prism phenomena. As an example, we consider a 2D 
PhC made of 5 lamellar grating layers (ei = 9, £2 = £ext = 1, di/d = 0.77, h = d, hi/d = 1/4). 
The isofrequency Bloch diagram of the structure is given in fig. 6 for X/d = 2. There is a zone 
of contra-propagating Bloch modes around ao = 1-6 {6 ~ 40°). The structure is illuminated 
by a monochromatic gaussian beam {w = lOA, = l.Q,X/d = 2). The map of the field 
is given in fig. 7, where it is seen that the shift of the transmitted beam is negative, the 
amplitude of the field on the upper and lower faces are given in fig. 8. The shift obtained by 
the direct numerical computation is A ~ —9.3, the shift obtained from the Bloch diagram 
is Af^/d = —3.3 whereas A^/c/ = —9.9. Once more, we find an excellent agreement between 
the direct computation and the effective theory, whereas the predictions of Bloch theory are 
quite false. 

Let us turn now toward a structure in which a guided contra-propagative mode do exist, 
this corresponds to the situation 2 in the proposition of the preceding section. In other 
words, there exists a mode with an horizontal wave vector. The parameters are the following: 
El = 7.84, di/d = 0.4:,hi/d = 0.7,h2/d = 0.3, X/d = 2.1, w = 50A, and there are 7 layers 
in the PhC. For this structure the contra-propagative mode is obtained for = 39°. Here, 



the parameter Cj which represents the behavior of (3 (a) ~ C-j^ — (x^ is obtained by a 
fitting of the results of the direct numerical computation, we obtained C7 = 0.164. We have 
plotted in fig. 9 [dj IS) versus dj [a — ao) where it is clearly seen that the shift converges 
towards a limit value. Using formula ( |22l) . we obtain A^/rf ~ —26 in fair agreement with 
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the numerical shift A/rf ~ —27. 

One should not think, however, that the non-renormalized Bloch diagram, i.e. that of 
the infinite periodic structure, cannot provide us with accurate results. It suffices to think 
of the homogenization regime, where the stack of gratings behaves as a stratified medium. 
For instance, we use a stack of 7 gratings {hi/d = 1/2, h = d,e = 2.1) and a beam with 
parameters: w = 25A, \/d = 4,00 = 40°. The field amplitudes on the upper and lower faces 
are given in fig. 10, where it can be seen that the oscillations are quite limited showing that 
we are indeed in the homogenization regime. The shift of the beam is A/c? = 8, and we have 
Ajj/d = 7.985 and A^/rf = 7.985. In that case, both predictions agree. This situation is 
due to the fact that the field inside the PhC can be represented by Bloch modes only. This 
situation may happen outside the homogenization regime. 



IV. CONCLUSION 



We have developped an effective medium approach to describe beam propagation inside 
a photonic crystal. This effective theory takes into account the evanescent waves, which 
are completely skipped if one uses only Bloch waves to describe wave propagation in the 
crystal. The importance of these evanescent waves are put into light by the computation 
of the shift of the transmitted beam. We show for some examples that the predictions 
obtained by using only the dispersion diagram may be false. These results emphasize the 
difference between the band theory for the Schrodinger equation (i.e. the propagation of 
electrons in periodic potentials), where the boundary of the crystal is irrelevant, and the 
scattering of electromagnetic waves by photonic crystals where boundary effects are of crucial 



importance. We have developed elsewhere 



14| some theoretical tools, that should hopefully 



permit to obtain a clearer insight into the role of evanescent waves. Work is in progress in 
that direction. 
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Appendix 1 



Lemma 

Let f be a real even function continuously differentiahle near ckq and such that f{ao) — 0. 
Then if f is square integrable near ao there holds 

= O (1) (23) 

Proof: 

Let us write: /(«) = / {a') + f (t) dt. Then (/ {a) - f {a')f < f (t) dtf and 



then by Cauchy-Schwarz inequahty: (/ (a) — / (a')) < \/a — a! \j (^J^, (/' (t)) dt^ and the 
proposition follows. 
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Appendix 2 



Applying the same reasoning as for the stratified medium, we get: 

Gt = -Nh J (a) dapN (a) da/Et 
Assuming that Uq has a first moment, we obtain, by Parseval equahty: 

J x\uo{x)fdx = \ J A (a) A (a)' da + J A^ (a) da^N {<^) da 
— J A^ (a) da^N (a) da 

and the angular shift is given by (cf. fig. 1): 

J A^ (a) daf^N (a) da 
1jB;Ii yj — — — 

JA^ (a) da 

The dispersion diagram is described locally by (3 — (f) (a) and, in the vicinity of the branch 
point ao, we can write from the lemma proved in Appendix 2: 

(a) — C ^ a^ — ccq, a>aQ 

the shift is then given by: 

jA^{a) . ^_ ^ da 3/4 . X 

tan^ = -C ^ °- — ^ -CJ^^J^^^^i^o^V [-\^kW (24) 

J A' {a) da A V4y ^ ' 

where ai = k sin 6i is a maximum of a — > /I (a)^ (a — aof^'^. 
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FIG. 1: sketch of the photonic crystal 
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FIG. 2: Basic cell of the photonic crystal used in the numerical experiments. 
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FIG. 3: Incident (a) transmitted (b) and reflected (c) fleld intensities for an s-polarized incident 
field on a ID structure. 
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(b) 
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FIG. 4: Incident (a) transmitted (b) and reflected (c) fleld intensities for a p-polarized incident 
field on a ID structure 
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(a) (b) (c) 

FIG. 5: Incident (a) transmitted (b) and reflected (c) field intensities for an s-polarized beam on 
a 2D PhC. 
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FIG. 6: Bloch diagram for the inverted contrast photonic crystal. 
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FIG. 7: Bloch diagram for the photonic crystal allowing contra-propagating modes. 
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x/d 



FIG. 8: Map of the electric field for the contra-propagating mode. 
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FIG. 9: Incident (a) transmitted (b) and reflected (c) field intensities for an s-polarized beam on 
a 2D PhC. 
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FIG. 10: Evolution of the shift with respect to the Fourier variable. 
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FIG. 11: Incident (a) transmitted (b) and reflected (c) field intensities in the homogenization 
regime. 
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